---
title: "WHAP 2019 Data Re-transform"
output:
  html_document:
    df_print: paged
---

The WHAP 2019 field data CSV on ServCat had significant errors when checked in
May 2023. These errors were corrected in analysis but the file with those
corrections has a different format than the current WHAP scripts use. This
notebook documents the process of correcting that.


```{r setup}
# Set up user directories
results_dir <- "results"
code_dir <- "src"
data_dir <- "data"

# Set up file names
corrected_csv_file_name <- "qdt2019c.csv"
servcat_csv_file_name <- "WHAP_2019_Quadrats_20211019.csv"
output_csv_file_name <- "WHAP_2019_Quadrats_20230525.csv"

required_packages <- c(
  "tidyverse", 
  "here"
)

# (optional) Initialize here::here()
here::i_am("src/whap_2019_data_retransform.Rmd")

# Load our (optional) function to make package loading more graceful
public_download_api <- 
  "https://ecos.fws.gov/ServCatServices/servcat/v4/rest/DownloadFile/"
try(source(paste0(public_download_api, "228091")))

if (exists("load_and_install_packages")) {
  load_and_install_packages(required_packages)
} else {
  message("Loading required packages: ", toString(required_packages))
  loaded_packages <- sapply(required_packages, 
         \(x) library(x, character.only = TRUE, quietly = TRUE)
  )
}

# (Optional) function for nicer messages in HTML RMarkdown docs
try(source(paste0(public_download_api, "227695")))
if (exists("make_errors_prettier") && !interactive()) make_errors_prettier()

# Make sure desired directories exist
if (!dir.exists(here(data_dir))) dir.create(here(data_dir))
if (!dir.exists(here(results_dir))) dir.create(here(results_dir))
```

## Load the input csvs

The corrected 2019 data was received from Emilio Laca via email on 10 May 2023.

The uncorrected 2019 data was pulled from ServCat on 26 April 2023.

```{r load-input-csv}
cli::cli_alert_info("Loading {corrected_csv_file_name}")
corrected_data <- read_csv(here(data_dir, corrected_csv_file_name))

cli::cli_alert_info("Loading {servcat_csv_file_name}")
servcat_data <- read_csv(here(data_dir, servcat_csv_file_name))

```

## Get event-level good data from ServCat data

Some of the columns on the ServCat data are already how we want them - the
issues are mostly with the taxon columns and measurements.

```{r extract-quadrat-data-from-SC}
quadrats_from_servcat <- servcat_data %>% 
  dplyr::select(GlobalID, LIT, unitName, subunitName, recordedBy, eventDate, 
                decimalLatitude, decimalLongitude, managementAction, 
                quadratSize, plantHeight, nSeedHeads, eventRemarks
  ) %>% 
  dplyr::distinct() %>% 
  dplyr::mutate(
    # The latitude and longitude fields are switched, just fixing them
    tempLatitude = decimalLongitude,
    decimalLongitude = decimalLatitude, 
    decimalLatitude = tempLatitude,
    tempLatitude = NULL
    ) 

head(quadrats_from_servcat)
```
## Get seedhead-level good data from corrected data

We only need a few columns from the seedhead data. First we will get the
taxonomy and stratum columns.

```{r extract-taxonomy-data-from-corr}
corrected_taxonomy <- corrected_data %>% 
  dplyr::select(GlobalID = OBJECTID, vernacularName = CommonName,
                taxonID = ITIS_TSN, stratum = Stratum
  ) %>% 
  dplyr::mutate(vernacularName = if_else(
    vernacularName == "Timothy",
    "Swamp Timothy",
    vernacularName
    ),
     stratum = if_else(
      stratum == "Med",
    "Medium",
    stratum
    ),
    stratum = stringr::str_to_lower(stratum)
  ) %>% 
  dplyr::distinct()

head(corrected_taxonomy)
```
### Transform seedhead-level measure data from corrected data

There are ten measure columns we want to pivot and get data from.

```{r pivot-seedhead-data-from-corr}
seed_head_cols <- names(corrected_data) %>% 
  grep("^seed_head", ., value = TRUE )

pivoted_measures <- corrected_data %>% 
  dplyr::select(GlobalID = OBJECTID, all_of(seed_head_cols)) %>% 
  tidyr::pivot_longer(
    seed_head_cols,
    names_to = "measurementType",
    values_to = "measurementValue",
    values_drop_na = TRUE
  ) %>% 
  dplyr::mutate(
    # Generate a GlobalID for the seed head measures
    GlobalID_seed = paste(
      GlobalID, 
      stringr::str_extract(
        measurementType,
        "(?!^seed_head)[1-5]" 
      ),
      sep = "_"
    ),
    # Replace measurement types with the names used in current data
    measurementType = if_else(
      stringr::str_detect(measurementType, "F2TLength"),
      "f2t_length_mm",
      "sh_length_mm"
    )
  ) 

head(pivoted_measures)
```

## Join all the data back together

We now have three datasets to stick together:

`quadrats_from_servcat`

`corrected_taxonomy`

`pivoted_measures`

```{r join-all-data}
full_corrected_data <- quadrats_from_servcat %>% 
  dplyr::inner_join(corrected_taxonomy, by = join_by(GlobalID)) %>% 
  dplyr::inner_join(pivoted_measures, by = join_by(GlobalID)) %>%
  # Just for convenience, put the columns in the original ServCat csv order
  dplyr::relocate(names(servcat_data)) %>% 
  dplyr::arrange(GlobalID, GlobalID_seed)

head(full_corrected_data)

```

## Checks on data completeness

I used GlobalID_seed to find records that appeared in the ServCat data but not
the full corrected data, and vice versa.  

These records appear in the full corrected data data but not the ServCat data. 
In  all cases I checked, it was because measurements were incorrectly labeled as 
F2TLength instead of regular seed head length (which is not measured in swamp
timothy).

`r setdiff(full_corrected_data$GlobalID_seed, servcat_data$GlobalID_seed)`

These records appear in the ServCat data but not the full corrected data. Mostly
these are records where the measurements are impossible, uncorrectable, or 
could not be resolved to a correct unit of measurement (like seed heads that 
were measured at 0.05 mm long).

`r setdiff(servcat_data$GlobalID_seed, full_corrected_data$GlobalID_seed)`

```{r write-out-to-csv}
write_csv(full_corrected_data, here(results_dir, output_csv_file_name))
```

